THE PROPAGATION AND ERUPTION OF RELATIVISTIC 
JETS FROM THE STELLAR PROGENITORS OF GRBs 



Weiqun Zhang and S. E. Woosley 

Department of Astronomy and Astrophysics^ University of California^ Santa Cruz, CA 

95064 

zhang@ucolick. org, woosley@ucolick. org 

and 
A. Heger 

Theoretical Astrophysics Group, T-6, MS B227, Los Alamos National Laboratory, Los 

Alamos, NM 87545 
Enrico Fermi Institute, University of Chicago, Chicago, IL 60637 

l(§2sn. org 
ABSTRACT 

New two- and three-dimensional calculations are presented of relativistic jet 
propagation and break out in massive Wolf-Rayet stars. Such jets are thought 
responsible for gamma-ray bursts. As it erupts, the highly relativistic jet core 
(3 to 5 degrees; F > 100) is surrounded by a cocoon of less energetic, but still 
moderately relativistic ejecta (F ^ 15) that expands and becomes visible at 
larger polar angles 10 degrees). These less energetic ejecta may be the origin 
of X-ray flashes and other high-energy transients which will be visible to a larger 
fraction of the sky, albeit to a shorter distance than common gamma-ray bursts. 
Jet stability is also examined in three-dimensional calculations. If the jet changes 
angle by more than three degrees in several seconds, it will dissipate, producing 
a broad beam with inadequate Lorentz factor to make a common gamma-ray 
burst. This may be an alternate way to make X-ray flashes. 

Subject headings: gamma rays: bursts — hydrodynamics — methods: numerical 
— relativity 
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1. INTRODUCTION 

It is now generally acknowledged that "long-soft" gamma-ray bursts (henceforth, just 
"GRBs") are a phenomenon associated with the deaths of massive stars. In addition to 
the observed association with star- forming regions in galaxies (Vreeswijk et al. 2001; Bloom, 
Kulkarni, & Djorgovski 2002; Gorosabel et al. 2003); "bumps" observed in the afterglows of 
many GRBs (Reichert 1999; Galama et al. 2000; Bloom et al. 2002; Garnavich et al. 2003; 
Bloom 2003); spectral features like a WR-star in the afterglow of GRB 021004 (Mirabal et 
al. 2002); and the association of GRB 980425 with SN 1998bw (Galama et al. 1998), there is 
now incontrovertible evidence that at least one GRB (030329) was accompanied by a bright, 
energetic supernova of Type Ic, SN 2003dh (Price et al. 2003; Hjorth et al. 2003; Stanek 
et al. 2003). Thus some, if not all GRBs, are produced when the iron core of a massive 
star collapses either to a black hole (Woosley 1993; MacFadyen & Woosley 1999) or a very 
rapidly rotating highly magnetic neutron star (Wheeler et al. 2000), producing a relativistic 
jet. Other alternatives in which the GRB occurs after the supernova (Vietri & Stella 1998) 
cannot simultaneously explain GRB 030329 and SN 2003dh. 

At the same time, the general class of high energy transients once generically called 
"gamma-ray bursts" has been diversifying. For some time, soft gamma-ray repeaters (SGRs) 
and "short hard gamma-ray bursts" have enjoyed a distinct status and, presumably, a sep- 
arate origin. In addition, we now have cosmological "X-ray flashes" ("XRFs"; Heise et al 
2001; Kippen et al. 2003), "long, faint gamma-ray bursts" (in't Zand et al. 2003), and lower 
energy events like GRB980425 (Galama et al. 1998). Is a different model required for each 
new phenomenon, or is some unifled model at work, as in active galactic nuclei, a model 
whose observable properties vary with its environment, the angle at which it is viewed, and 
perhaps its redshift? 

The answer is probably "both" . Even within the conflnes of the coUapsar model, there 
are Types I, II and III (Heger et al. 2003). Type I happens only in massive stars that make 
their black holes promptly (MacFadyen & Woosley 1999). This is what people generally have 
in mind when they use the word "coUapsar" . In Type II though, a similar mass black hole 
forms by fall back and the burst lasts much longer (MacFadyen, Woosley, & Heger 2001). 
Type III, happens only at very low metallicity and requires Pop III stars as its progenitors 
(Fryer, Woosley, & Heger 2001). The event is very energetic, but highly redshifted. These 
other models may be particularly appropriate for unusually long bursts. A jet that wavered in 
its orientation in any of these events would emerge with a greater load of baryons. So would 
a jet whose power declined in a time less than that required to traverse the star. A jet in a 
blue supergiant would become choked because the power of the central engine has decreased 
substantially when the jet is still deep inside the star. Thus any sort of coUapsar happening 
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in a blue supergiant instead of a Wolf-Rayet star would produce a different sort of transient. 
Finally, it may just be that jets are just born with different baryon loadings, perhaps relating 
to the mass of the presupernova star and its distribution of angular momentum and magnetic 
field. Our present understanding of nature allows for all these solutions. 

However, it also seems to reason that not all jets are the same, even for Type I coUapsars, 
and that even if they were, different phenomena would be seen at different angles. We are 
thus motivated to consider the observational consequences of highly relativistic jets as they 
propagate through, and emerge from massive stars. What would they look like if seen from 
different angles? In particular, what is the distribution with polar angle of the energy and 
Lorentz factor? Might softer, less energetic phenomena be observed off axis and with what 
frequency? 

Jets inside massive stars have been studied numerically in both Newtonian (MacFadyen 
& Woosley 1999; MacFadyen, Woosley, & Heger 2001; Khokhlov et al. 1999) and relativistic 
simulations (Aloy et al. 2000; Zhang, Woosley, & MacFadyen 2003, henceforth Paper 1) 
and it has been shown that the coUapsar model is able to explain many of the observed 
characteristics of GRBs. These previous studies have also raised issues which require further 
examination, especially with higher resolution. For instance, the emergence of the jet and its 
interaction with the material at the stellar surface and the stellar wind could lead to some 
sort of "precursor" activity. The cocoon of the jet will also have different properties from 
the jet itself and shocks within the cocoon or with external matter could lead to 7-ray and 
hard X-ray transients (Ramirez- Ruiz, Celotti, & Rees 2002). There is also the overarching 
question of whether jets, calculated in two dimensions with assumed axial symmetry, are 
stable when studied in three dimensions. 

Here we examine in two- and three-dimensional numerical studies, the interaction of 
relativistic jets with the outer layers of the Wolf-Rayet stars thought responsible for GRBs. 
We follow the emergence of the jets for a sufficient length of time to ascertain the properties 
of the cocoon explosion that surrounds the GRB and we study the dependence of the results 
on the dimensionality of the calculation. Such explosions are visible to a much greater angle 
and to a much larger number of observers. The observed phenomenon could be a hard X-ray 
flash or a weak GRB. We also find that the properties of the emergent jet are quite sensitive 
to whether the jet maintains its orientation (to within a critical angle) for durations of 10 s or 
so, the time it takes the jet to traverse the star. Jets that waiver may make hard transients 
of some sort, but not GRBs. 
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2. Stellar Model and Numerical Methods 

2.1. Progenitor Star 

A coUapsar is formed when the iron core of a rotating massive star coUapses to a black 
hole and an accretion disk. The interaction of this disk with the hole, through processes 
that are still poorly understood, produces jets with a high energy to mass ratio. Here we are 
concerned primarily with the propagation of these relativistic jets, their interactions with 
the star and the stellar wind, and the observational implications, and not so much with how 
they are born. In fact, our results here would be the same if the jet were produced by some 
other means (Wheeler et al. 2000) inside the same star. 

Unlike in our previous study (Paper 1), the initial stellar model is taken directly from 
a presupernova star with a very finely zoned surface (Table 1; Fig. 1). Because we are 
only interested in events that happen outside about 10^° cm on a time scale of, at most, 
tens of seconds, modification of the initial star by core collapse and rotation is negligible. 
Specifically, our initial model is a 15 helium star calculated by Heger & Woosley (2003). 
This helium star has been evolved to iron-core collapse while following the transport of 
angular momentum and including the effects of rotational mixing as in Heger, Langer, & 
Woosley (2000). The initial star was a nearly pure helium star rotating rigidly with a surface 
velocity equivalent to 10% Keplerian. This is a typical value in the study of Heger, Langer, 
& Woosley (2000). Presumably the star lost its envelope to a companion early in helium 
burning. Further mass loss and the transport of angular momentum by magnetic fields were 
ignored so as to form a big iron core that would very likely collapse to a black hole and 
give ample angular momentum for a disk. The radius of the helium star at onset of collapse 
is 8.8 X lO"*^^ cm and the surface of the star is very finely zoned (surface zoning < lO^-'^g). 
Despite the fact that angular momentum transport was followed in the initial model to show 
its suitability as a coUapsar (Fig. 1), rotation plays no role in the present study. 

Outside the star, the background density, which comes from the stellar wind, is assumed 
to scale as R~^. The actual value of this density, unless it is very high, plays no role in 
determining the final answer at the small scales considered here (< 2 x 10^^ cm), but finite 
value is needed in order to stabilize the code. We assumed a value at i? = 10^^ cm of 
5 X 10~^^gcm~^, where R is the distance from the center of the star. This corresponds to 
a mass loss rate of ~ 1 x lO~^M0yr~^ for a wind velocity of ^ 1000 km s~^ at 10^^ cm. 
This mass loss rate is rather low for typical Wolf-Rayet stars of this mass in our galaxy 
(Nugis & Lamers 2000), but might be more appropriate for the low metallicity in the GRB 
neighborhood (Crowther et al. 2002). A small value was also taken in order to be consistent 
with the assumed progenitor structure that was calculated with zero mass loss. 
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2.2. Computer Code 

We employed a multi-dimensional relativistic hydrodynamics code that has been used 
previously to study relativistic jets in the coUapsar environment (Paper 1). Briefly, our code 
employs an explicit Eulerian Godunov-type shock-capturing method (Aloy et al. 1999). The 
governing equations of relativistic hydrodynamics with a causal equation of state can be 
written as a system of conservation laws for rest mass, momentum, and energy. To solve 
the equations numerically, each spatial dimension is discretized into cells. Using the method 
of lines, the multi-dimensional time-dependent equations can then be evolved by solving 
the fluxes at each cell interface. In order to achieve high order accuracy in time, the time 
integration is done using a high order Runge-Kutta scheme (Shu & Osher 1988), which 
involves prediction and correction. An approximate Riemann solver (Aloy et al. 1999) using 
the Marquina's algorithm is used to compute the numerical fluxes from physical variables: 
pressure, rest mass density, and velocity at the cell interface. The values of the physical fluid 
variables at the cell interface are interpolated using reconstruction schemes (e.g., piecewise 
parabolic method, Colella & Woodward 1984). This reconstruction procedure ensures high 
order accuracy in spatial dimensions. The conserved variables: rest mass, momentum and 
energy are evolved directly by the scheme. In each time step, physical variables, such as 
pressure, rest mass density and velocity, which are necessary in calculating the numerical 
fluxes, can be recovered from conserved variables by Newton-Raphson iteration. The code 
can operate in Cartesian, cylindrical, or spherical coordinates. Approximate Newtonian 
gravity is implemented by including source terms in the equations. Total energy, which 
includes kinetic and internal energy, in the laboratory frame are not exactly conserved due to 
gravitational potential energy and the way of implementing gravity as source terms. However, 
the code conserves energy to machine precision when gravity is turned off. Furthermore, the 
gravitational potential energy is negligible for relativistic material in our calculations. For 
instance, assuming a 10 Mq central point mass, a fluid element with a velocity of 0.5 c at the 
inner boundary of our computational grids, 10^° cm (§ 3.1, § 4.1), has potential energy of less 
than 0.1% of its kinetic energy. Thus we expect that our calculations conserve energy to an 
accuracy of better than 0.1%. In simulations present in this paper, A gamma-law equation 
of state with 7 = 4/3 is used for the simulations presented in this paper. 
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3. Two-Dimensional Models 

3.1. Model Set Up and Definitions 

The mass interior to 1.0 x 10^^ cm is removed from the presupernova star and replaced by 
a point mass. No self gravity is included. This should suffice since we are studying phenomena 
that happen on a relativistic time scale and the speed of sound is very sub-luminal. While 
jets presumably go out both axes, we follow here only one, assuming symmetry in the other 
hemisphere. Jets are injected along the rotation axis (the center of the cylindrical axis of 
the grid) through the inner boundary. Each jet is defined by its power (excluding rest mass 
energy), its initial Lorentz factor, Fq, and the ratio of its total energy (excluding rest 
mass energy) to its kinetic energy, /q. Our previous studies in Paper 1 have shown that a jet 
that starting with a half-opening angle of 20° and a high Lorentz factor, F ^ 50 at 2000 km, 
will be shocked deep inside of the star. By the time it reaches 10^^ cm, a jet should have a 
large ratio of internal energy to rest mass, a half-opening angle of about 5°, and a Lorentz 
factor, F ^ 5 — 10. 

Though GRBs observationally have highly variable properties, the jet power was taken 
here to be constant for the first 20 s, then turned down linearly during the next 10 seconds. 
That is, during the interval 20 to 30 s, the power scaled as (30 s — t)/10s decreasing to zero 
at 30 s. During the declining phase, the pressure and density remained constant, and the 
Lorentz factor was calculated from the internal energy, density, and power. After 30 seconds, 
a pure outfiow boundary condition was used for the inner boundary. The axis of the jet is 
defined as the z-axis in all cases and the jets were initiated parallel to that axis in a region 
that subtended a half-angle of 5 degrees as viewed from the origin. 

Four models were calculated which span a range of energies and Lorentz factor (Table 2). 
In Model 2A, a total energy of ^ 5 x 10^^ erg is injected, comparable to the results of Frail 
et al. (2001), but larger than the results of Panaitescu & Kumar (2001) by an order of 
magnitude (Table 3). In Models 2B and 2C, higher and lower energy deposition rates were 
employed. For the parameters chosen, a jet with an initial Lorentz factor of 5 or 10 at 
10^° cm should have a final Lorentz factor of ~ 180 — 200 if all internal energy is converted 
into kinetic energy. The zoning employed in the models is given in Table 4. Typically over 
three million zones were used. 

Model 2C used an extended grid in the r-direction, which allowed the greater lateral 
expansion of the jet to be followed at late times. This was necessary because of its lower 
energy and greater expansion during the time it took to traverse the extent of the z-grid. 
Model 2T used conditions like those of Model 2B, but a grid that was both smaller and 
coarser, chosen to be equivalent to that used in the three-dimensional calculations described 
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later in the paper (§ 4). 

3.2. Results in Two Dimensions 

The relativistic jet begins to propagate along the z-axis shortly after its initiation. 
In agreement with previous studies (Aloy et al. 2000, and Paper 1), the jet consists of a 
supersonic beam, a shocked cocoon, a bow shock, and is narrowly collimated. Some snapshots 
of Models 2A, 2B, and 2C are given in Figs. 2, 3, and 4. The resulting "equivalent isotropic 
energies" for matter exceeding a certain Lorentz factor are also given, at 70s, for the same 
three models in Figs. 5, 6, and 7. The latter figures also give the estimated terminal Lorentz 
factor assuming that internal energy along the radial line of sight is converted into kinetic 
energy. The equivalent isotropic energy is defined as that energy an isotropic explosion 
would need in order to give the calculated fiux of energy along the line of sight. It is a way 
of mapping the highly asymmetric two-dimensional results into equivalent one-dimensional 
models. The change in Lorentz factor between 70s and infinity is not very great, except 
in situations where matter that was highly relativistic to begin with receives an additional 
boost in its frame by expansion. The duration of the calculation was set by how long it took 
the jet to reach the end of the simulation grid and would be costly to increase, but as the 
figures show, especially in Model 2C, there was still an interesting amount of internal energy 
at 70s. 

The fractions of energies inside a certain angle to the total energies on the grid are 
given in Fig. 8. In all cases studied the high Lorentz factor characteristic of common GRBs 
is confined to a narrow angle of about 3 to 5 degrees with a maximum equivalent isotropic 
energy in highly relativistic matter along the axis of ^ 3 x 10^^ to 3 x 10^^ erg. At larger 
angles there is significant energy, though, and Lorentz factor, F ^ 10 to 20. At an angle of 
10 degrees for example, the equivalent isotropic energy in matter with F > 20 is ^ 10^^ erg 
in Model 2A and even 10^^ erg in Model 2C. At larger angles there is less energy, but still the 
possibility of low power transients of hard radiation. Model 2B did not eject much material 
with F > 10 at angles ^ 10 degrees. The equivalent isotropic energy at larger angles (> 2°) 
for all three models can be fit well by a simple power-law, Eq x (^/2°)~^erg, where £"0 is 
1.5 X 10^^, 4.5 X 10^^ and 6.8 x 10^^ for models 2A, 2B and 2C, respectively. The ratios 
of the values of £"0, 1.5 : 4.5 : 0.68, are very close to those of the energy deposition rates, 
1.0 : 3.0 : 0.5. Inside 2°, the distributions of energy and Lorentz factor are roughly fiat, 
^ Eq. The values of Eq can be roughly estimated from the total injected energy and the 
percentage of the energy in the jet core. For example, about 40% of the total energy on the 
grid is contained inside 2 degrees for Model 2A (Fig. 8). About 5 x 10^-^ erg (Table 2) was 
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injected for both jets in Model 2A. One can estimate that the equivalent isotropic energy 
inside 2 degrees is ^ 1.6 x 10^^ erg for Model 2A. More models need to be calculated to find 
how the distributions of energy and Lorentz factor at breakout depend on initial conditions. 

Table 3 gives the energies for various components of the ejecta at a time 40 s after the 
jet was initiated. There still remains considerable internal energy available for conversion to 
kinetic at this time (which was chosen as a time when most of the relativistic matter was still 
on the grid), but this will affect mostly the matter with high Lorentz factor, F > 10. The 
results show that most of the energy injected as jets at 10^° cm emerges in matter that is still 
relativistic. Only a small fraction goes into sub-relativistic expansion, nominally v < 0.5 c 
and consequently, only a little of the jet energy is used to blow up the star. If there were no 
other energy sources the kinetic energy of the supernova would be < 2 x 10^^ erg. 

The total energy of relativistic ejecta is also useful for comparison to radio observations 
of the afterglows of GRBs. For example, Li & Chevalier (1999) have placed limits on the 
total energy of ejecta with v > 0.5 c in SN 1998bw. The limit, 3 x 10^^ erg, may be an 
approximate estimate of the actual energy. The corresponding energy for Model 2C here is 
four times larger (Table 3). It seems likely that lower energy models than Model 2C could 
be constructed that would still give high Lorentz factors and equivalent isotropic energies 
in a narrow range of angles around the polar axis. That is, GRB 980425 may have been a 
harder, more energetic GRB seen off axis. However, the low energy of that burst, some five 
orders of magnitude less that that expected for a centrally observed GRB (Fig. 4) remains 
surprising. Either GRB 980425 was observed at a polar angle larger than 15 degrees, or the 
burst was weaker at all angles because of baryon loading (§ 1; Woosley, Eastman, & Schmidt 
1999). 

The possibility that the off- axis emission from material with F ~ 10 to 20 corresponds 
to XRFs is discussed in § 5. 

In 2D simulations with cylindrical coordinates, there is an imposed symmetry axis of the 
coordinate system. It is very important to repeat these calculations in 3D to ensure that our 
2D results are valid and to examine three-dimensional jet instabilities. Three-dimensional 
calculations are, however, computationally expensive. So we have to use lower resolution 
for 3D calculations. In order to compare 2D and 3D results with the same resolution and 
identical initial parameters, we did a lower resolution 2D run. Model 2T (Tables 2 and 4). A 
comparison of two- and three-dimensional results will be discussed in § 4. It is also interesting 
to compare Model 2B and 2T, which had identical parameters but in different resolution. 
Qualitatively the results are similar. In both models, the jet emerges from the star with a 
cocoon surrounding the jet beam and a dense "plug" at the head of the jet (Fig. 9). There 
are, however, noticeable differences, as we expected. In particular, it takes less time for the 
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jet in Model 2T to emerge from the star than that in Model 2B, because they have different 
numerical viscosity. The larger numerical viscosity in Model 2T is due to its lower resolution. 
Despite morphologic differences between Model 2B and 2T (Fig. 9), the energetics of these 
two models are quite similar (Table 5, Fig. 10). The distributions of equivalent isotropic 
energy versus angle for the jet core (^ 3 degrees) in Models 2B and 2T are very similar. 
However, Model 2B has more mildly relativistic material (2 < F < 10) in its cocoon than 
Model 2T. 



4. Three-Dimensional Models 

For the three-dimensional models, the same helium star was remapped into a three- 
dimensional Cartesian grid. Because of the greater cost of computations, these studies 
covered only an interval of 10 s, adequate to watch the jet propagate, develop a cocoon and 
break out, but not long enough to study the the expansion after break out to any great 
extent. 



4.1. Model Definitions 

The parameters of the jet, its power, Lorentz factor and energy loading, were all identical 
to those of Model 2B, that is F = 5, ^ = 3 x 10^° erg s"^ and fo = 40 (Table 6). The grid 
employed in all cases was Cartesian with 256 zones each along the x- and y-axes and 512 
along the z-axis (jet axis) (Table 7). 

The initial model and the initiation of the jet in Model 3 A were perfectly symmetric 
with regard to the axis of the jet. The perfect symmetry was maintained in the calculation 
because our numerical scheme did not break any symmetry. In order to break the perfect 
symmetry of the cylindrical initial conditions, an asymmetric jet was injected in Models 3BS 
and 3BL. At the base of the jet in Model 3BS, pressure and density were 1% more than those 
of Model 3A if ^ > tan ax, and 1% less otherwise, where a = 40 degrees. Hence the jet in 
Model 3BS has a ±1% imbalance in power. In Model 3BL, a ±10% imbalance in power was 
employed. 

Several other models were calculated to explore the collimation properties of non-radial 
jets that precessed. These jets were initiated at the same 10^° cm inner boundary with a 
half-angle of 5 degrees as measured from the center. However, they were given a non-z 
(symmetry axis) component of momentum that would have resulted, in a vacuum, in a 
propagation vector inclined to the z-axis by 3 degrees (Model 3P3), 5 degrees (Model 3P5), 
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and 10 degrees (Model 3P10). The jet precessed around the z-axis with a period, in aU cases 
of 2 s. This period was chosen to be short compared with the time it took the jet to traverse 
the remaining star between 10^^ cm and the surface, yet long enough that the jet would still 
be distinct from a cone. 



4.2. Breakout in Three Dimensions 

As expected. Model 3A closely resembles Model 2T (Fig. 11). The initial stellar model, 
the zoning (Tables 4 and 7), and the jet parameters (Tables 2 and 6) were all the same. In 
particular, the initial conditions for the jet have cylindrical symmetry in both cases. It is 
gratifying that the answer is insensitive to the dimensionality of the grid, which in Model 2T, 
was two-dimensional cylindrical and, in Model 3A, was three-dimensional Cartesian. This 
implies that the bulk properties of jets calculated in Paper 1 - including collimation and 
modulation - would be the same in 3D. 

However, Model 3 A, by assuming a jet that initially has perfect cylindrical symmetry 
does not fully exercise the 3D code. Aside from numerical noise, the equations conserve 
the 2D symmetry of the initial conditions, even on a 3D grid. Fig. 11 also compares, at 
breakout, the properties of jets that were, initially, nearly identical. In particular. Model 
3BS differed from 3 A only in a 1% asymmetry in input energy from one side of the jet to 
the other, yet the structure of the emergent jet and cocoon is strikingly different. In Model 
3BL with a 10% asymmetry the difference is even more striking. The plots show a cross 
section of the Cartesian grid along the initial jet axis in the x = plane. Because of the 
40 degree offset (§ 4.1), there is somewhat more energy in the top half of the figure than 
the bottom. Though both jets retain high Lorentz factors in their cores, similar to Models 
2T and 2B, the cocoon explosion is a little earlier and larger on the top. More dramatic 
is the difference in the high density "plug" among Models 3A, 3BS, and 3BL. In the latter 
two where the 2D symmetry was mildly broken, the plug has a much lower density and is 
not prominent. In 2D models and symmetric 3D Model 3 A, the plug is held by a concave 
surface of the highly relativist ic jet core. Because of the imposed axisymmetry, the plug 
cannot easily escape and is pushed forward by the jet beam. Whereas the story of the plug 
is different in asymmetric 3D Models 3BS and 3BL. In these models, where the 2D symmetry 
was initially mildly broken, instabilities will develop. The forward bow shock is no longer 
symmetric and even the stellar material on the axis is pushed sideways. More importantly, 
the head of the highly relativist ic jet beam is also asymmetric and does not have a concave 
surface to "hold" the plug. A movie of these runs shows the plug forming, slipping off to 
the side, then forming again. The plug in these models has a much lower density and is not 
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prominent. The presence or absence of this plug may have important imphcation for the 
production of short-hard gamma-ray bursts in the massive star models (Paper 1; Waxman 
& Meszaros 2003). 

4.3. Stability of the Jet 

Given the ability to model jets in 3D, we undertook a study to test their survivability 
against non-radial instabilities. Three jets were introduced on the standard grid (Table 7) 
with inclinations to the radial of 3, 5, and 10 degrees (Table 6). These jets were made to 
precess with a period equal to 2 s, i.e., a non-trivial fraction of the time it takes the jet to 
traverse the star and break out. The results are shown in Fig. 12. 

With an angle of 3 degrees, the jet in Model 3P3 escapes the star with its relativistic 
flow at least partly intact. Though Fig. 12 shows Lorentz factors of only F ^ 10 — 20, the 
internal energy is still very large and some of these regions will attain F ~ 180 if they expand 
freely. However, there is some "baryon-poisoning" . No clear line of sight exists along a radial 
line for the most energetic material even in the jet beam. Perhaps a jet with intermediate 
Lorentz factor would emerge. We have not yet followed these calculations to large radius 
(as in Figs. 2, 3, and 4) because the 3D grid was not large enough. Larger calculations are 
planned. 

At larger angles Model 3P5 and especially 3P10 show the break up of the jet. Because 
there is no well focused highly relativistic jet beam, especially in Models 3P5 and 3P10, 
more baryon mass is mixed into the jet. The Lorentz factor of the jet will decrease due to 
dissipation. And it will be very difficult for these jets to make a common GRB. Again some 
sort of hard transient might be expected, especially after all the internal energy converts and 
F rises, but for Model 3P10 there will be no GRB of the common variety. In our simulations, 
we find that the critical angle for jet precession is about 3 degrees. One would expect that 
the constraint on the angle of precession will be reduced if the jet bears more power or is 
powered longer. 

Is two seconds a reasonable period for the jet to precess? There are many uncertainties, 
but the gravitomagnetic precession period of a black hole surrounded by an accretion disk 
is estimated by Hartle, Thorne, & Price (1986) to be 

Tgm = 9OO(Mdisk/O.OOlM0)-i(Mbh/3M0)-i/2(^/3QQj^^)2.5g^ 

where M^isk is the mass of the accretion disk, Mbh is the mass of the hole and r is the 
radius of the disk, r = 300 km is a reasonable value for the radius of the disk, which 
depends on the distribution of angular momentum of the star (MacFadyen & Woosley 1999). 
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The dependence of the mass in the disk on the alpha-viscosity has been explored in both 
numerical simulations (MacFadyen & Woosley 1999) and analytic calculations (Popham, 
Woosley & Fryer 1999). For currently favored value of the alpha- viscosity, ~ 0.1, the disk 
mass is about 0.001 Mq. For the disk mass to approach 1 the viscosity would need to be 
a < 0.001. Given the low mass of the accretion disk and the radius of the disk, ^ 300 km, 
precession periods as small as two seconds are unlikely, but the black hole may be kicked 
or instabilities deeper in the region of the star not modeled here could give the jet some 
non-radial component. If so. Fig. 12, shows an alternate way in which baryon-loaded jets 
could give rise to less relativist ic mass ejected and slower moving jets far away from the 
star(10 < r < 100). Softer transients like XRFs and sub-luminous GRBs like GRB 980425 
could result. 



5. Discussion 

Our calculations show, so long as the orientation of the jet does not waiver on time 
scales of order seconds by angles greater than 3 degrees, that a relativist ic jet can traverse 
a Wolf-Rayet star while retaining sufficient energy and Lorentz factor to make a GRB. This 
conclusion is robust in three dimensions as well as two. 

As it breaks out, the jet is surrounded by a cocoon of mildly-relativistic, energy-laden 
matter. As internal energy converts to kinetic, matter with 10^^ to 10^^ erg of equivalent 
isotropic energy moving with Lorentz factor F > 20 is ejected to angles about three times 
greater than the GRB. That is, whatever transient the cocoon gives rise to will be an order 
of magnitude more frequently observable in the universe, but two orders of magnitude less 
energetic than a GRB. Weaker transients can of course be obtained at still larger angles and 
there is considerable diversity in the models (Figs. 5, 6, and 7). 

The isotropic equivalent rest mass to 10^^ - 10^^ erg with F 20 is 3x 10"^ to 3x 10"^ M©. 
This matter would give up its energy upon encountering 1/F times its rest mass. For an 
assumed mass loss rate for the Wolf-Rayet progenitor of 10~^ yr~^ and speed, 1000 kms~^, 
the energy will be radiated at about 3 x 10^^ to 3 x 10^^ cm. Emission from this deceleration 
will have a duration ^ R/ (2F^c) ^ 10 - 100 s. Variations of a factor of 10 are easily achievable 
by varying the mass loss rate or Lorentz factor. 

Unlike the central jet, we do not see evidence for large scale variation in the Lorentz 
factor at lower latitude and so it would be hard to make the transient emission by internal 
shocks. We thus attribute the emission to external shocks and invoke a variety of Lorentz 
factors emitting within our light cone in order to explain the spread in wavelength. 
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Might there be observable counterparts to these large angle, low Lorentz factor ex- 
plosions? Typically, they have too low a Lorentz factor to make common GRBs (Guetta, 
Spada, & Waxman 2001). Tan, Matzner, & McKee (2001) present arguments however that, 
by external shock interaction with the progenitor wind, a hard transient of some sort should 
result. 

Our model predicts a correlation between Ei^o^ the Lorentz factor, and the angle between 
the polar axis and the observer (Fig. 5, 6 and 7). Roughly speaking, the GRB outflows 
have a narrow highly relativistic jet beam and a wide mildly relativistic jet wing. Recent 
observations and afterglow modeling support this non-uniform jet model (Berger et al. 2003; 
Zhang et al. 2003; Salmonson 2003; Rossi, Lazzati, & Rees 2002). If, as seems reasonable, the 
Lorentz factor is, in turn, correlated with the peak energy observed in the burst one expects 
a continuum of high energy transients spanning the range from X-ray afterglows (keV) , to 
hard X-ray transients (tens of keV) , to GRBs (hundreds of keV) . Observations reviewed by 
Amati et al. (2002) show such a correlation for bursts with E'peak from 80 keV to over 1 MeV 
and Lamb (private communication) flnds that the relation extends to the lower energies. 

In particular, XRFs form a new class of X-ray transients having a duration of order 
minutes and properties that in ways resemble GRBs (Heise et al 2001; Kippen et al. 2003). 
XRFs are probably many phenomena (Areflev, Priedhorsky, & Borozdin 2003) and it could 
be that, especially some of the longer ones, have alternate explanations (see introduction). 
However, we have felt for some time (Woosley, Eastman, & Schmidt 1999; Woosley & Mac- 
Fadyen 1999; Woosley 2000, 2001, Paper 1) that many XRFs are the oflF-axis emissions of 
GRBs, made in the lower energy wings of the principal jet. We have not calculated the 
expected spectra for any of our models, only Lorentz factors and energies. However, it is 
reasonable that matter moving with F ^ 20 would make a transient softer than a GRB and 
harder than a few keV. Larger amounts of energy are emitted at a larger range of angles for 
slower but still relativistic ejecta (F = 5, e.g.. Fig. 5). 

If our model is valid, XRFs and GRBs should be continuous classes of the same basic 
phenomenon sharing many properties. They should have a similar spatial distribution to 
GRBs because they are essentially the same sources. However, because they are much less 
luminous, their log N-log S distribution would not exhibit the same roll over attributed in 
GRBs to seeing the "edge" of the distribution (e.g., Fishman & Meegan 1995). Their median 
redshift should be considerably smaller, certainly less than 1. This is consistent with the 
low redshifts inferred for the host galaxies of two XRFs by Bloom et al. (2003). XRFs may, 
most frequently, be seen in isolation and will be characterized by softer spectra, but there 
would also be an underlying XRF in every GRB since the emission of the mildly relativistic 
cocoon material is beamed to a larger angle that includes the poles. In some cases these 
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XRFs might be seen as precursors or extended hard X-ray emission foUowing a common 
GRB. They should be associated with supernovae. Indeed XRFs may more frequently serve 
as guideposts to jet-powered supernovae than GRBs, especially the nearby ones. 

Though considerable variation is expected, our calculations (e.g., Fig. 5) suggest that 
XRFs are typically visible at angles about three times greater than GRBs and hence to ten 
times the solid angle. However, their energy, a few percent of GRBs, implies that a flux- 
limited sample could observe GRBs out to roughly ten times farther, implying that XRFs 
would be about 1% as frequent in the sample. The actual value is detector sensitive, but 
may be larger. 

One should also keep in mind the possibility that XRFs do not accompany GRBs, in any 
direction, but are a result of baryon loading of the central jet. Three-dimensional instabilities 
(Fig. 12) may play an important role in this and are being explored. 
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Table 1. Properties^ of the Progenitor Model 





m 


r 


J 




(Mo) 


(cm) 


(erg s) 


iron core 


1.95 


2.58x10^ 


1.23x10^° 


Si core 


2.61 


4.99x10^ 


2.35x10^° 


Ne/Mg/0 core 


2.95 


7.06 xlO^ 


2.68x10^° 


C/0 core 


8.56 


5.16x10^ 


2.24x10^^ 


star / He core 


15.00 


8.80x10^° 


1.00x10^2 



^Enclosed mass, m, radius, r, and enclosed angu- 
lar momentum, of the progenitor model at core 
collapse on the outer boundaries of the indicated 
cores. 
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Table 2. Parameters of 2D-Models 













rp^e 








Model 


(10^° ergs-1) 


J- 




(10^° cm) 


(s) 


(s) 


(IQSi erg) 


Grid'^ 


2A 


1.0 


10 


20 


1.0 


20 


10 


5.0 


Normal 


2B 


3.0 


5 


40 


1.0 


20 


10 


15.0 


Normal 


2C 


0.5 


5 


40 


1.0 


20 


10 


2.5 


Large 


2T 


3.0 


5 


40 


1.0 


20 


10 


15.0 


Coarse 



^Energy deposition rate per jet 
^Initial Lorentz factor 

^Initial ratio of the total energy (excluding the rest mass) to the kinetic energy 

^Low-z boundary, where the jet was injected 

^During this period, the jet power remained constant. 

^During this period, the jet power was turned down linearly. 

^Total energy injected during the calculation 

^See Table 4 for details 
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Table 3. Energetics of 2D-Models^ 



Model 


Energy input 


Total on grid'' 


w > 0.5c 


r > 2 


r > 10 




(IQSi erg) 


(IQSi erg) 


(IQSi erg) 


(10^1 erg) 


(IQSi erg) 


2A 


5.0 


4.66 


3.64 


3.56 


3.32 


2B 


15 


14.36 


12.46 


12.22 


11.16 


2C 


2.5 


1.93 


1.20 


1.13 


0.92 



^Energies evaluated for the entire star, both jets at 40 s. 

^The calculations conserve energy to an accuracy of better than 0.1%. How- 
ever, some energy has left the computational grid. 
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Table 4. Zoning in 2D-Models (10^° cm) 



Model 




r2 






Z2 


Z3 


r 


z 




An 


Ar2 


Ars 


Azi 






zones 


zones 


2A 


10.0 


20.0 


60.0 


10.0 


20.0 


200 


1500 


2275 




0.01 


0.04 


0.16 


0.01 


0.04 


0.16 






2B 


10.0 


20.0 


60.0 


10.0 


20.0 


200 


1500 


2275 




0.01 


0.04 


0.16 


0.01 


0.04 


0.16 






2C 


10.0 


20.0 


200.0 


10.0 


20.0 


200 


2375 


2275 




0.01 


0.04 


0.16 


0.01 


0.04 


0.16 






2T 


0.64 


3.84 




13.8 






256 


512 




0.01 


0.05 




0.025 











^The first row for each model gives tlie upper bound of distance (r, 
or z) for wliicli tlie zoning indicated in the second row was employed. 
All are measured in units of 10^° cm. r starts at 0; z starts at 1.0 x 
10^° cm. The zoning in Models 2A and 2B was identical. 
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Table 5. Energetics of 2D-Models 2B and 2T 



Model 


Energy input 


Total energy 


-y > 0.5c 


F > 2 


r > 10 




(10^1 erg) 


(10^^ erg) 


(10^^ erg) 


(10^1 erg) 


(10^^ erg) 


2B 


3 


2.46 


1.44 


1.31 


0.51 


2T 


2.85 


2.42 


1.50 


1.13 


0.52 



^Energies for Model 2T are evaluated for the entire computational grid at 
9.5 s, whereas, energies for Model 2B are evaluated for only part of the whole 
grid at 10 s. Note that the grid of Model 2B is larger than that of Model 2T 
(Table 4). We only consider part of the grid for Model 2B so that the same 
physical region in the two models are discussed. The bow shocks of the jets are 
at similar radius for the two models at the chosen moment (Fig. 9). Also note 
that we only consider one hemisphere. 
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Table 6. Parameters of 3D-Models 





E 








a a 
Up 


Period*^ 


Model 


(10^° ergs"^) 


To 


/o 


(10^° cm) 


(degrees) 


(s) 


3A 


3.0 


5 


40 


1.0 







3BL'= 


3.0 


5 


40 


1.0 







3BS^ 


3.0 


5 


40 


1.0 







3P3 


3.0 


5 


40 


1.0 


3 


2 


3P5 


3.0 


5 


40 


1.0 


5 


2 


3P10 


3.0 


5 


40 


1.0 


10 


2 



^ Angle between the axis of the injected jet and the z-axis 

^Period of precession in Models 3P3, 3P5, and 3P10 

^Models 3BL and 3BS were hke Model 3A (and Model 2B) but in- 
cluded an asymmetric energy input at the base of 10% in Model 3BL 
and 1% in Model 3BS. 



-24- 



Table 7. Zoning in 3D-Models (10^° cm) 





X2 


yi 


1/2 




X 


y 


z 


Axi 


AX2 


A|/i 


A?/2 


Azi 


zones 


zones 


zones 


0.64 


3.84 


0.64 


3.84 


13.8 


256 


256 


512 


0.01 


0.05 


0.01 


0.05 


0.025 









^The first row for each model gives tlie absolute value of 
the upper bound of distance (x, ^, or z) for which the zoning 
indicated in the second row was employed. All are measured in 
units of 10-^^ cm. \x\ and \y \ start at 0; z starts at 1.0 x 10-^^ cm. 
All three-dimensional studies used this same zoning. 
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Fig. 1. — Structure of the star at the onset of core coUapse as a function of enclosed mass, m. 
Panel A gives mass fractions of the dominant species, Panel B density (p) and temperature 
(T) stratification, and Panel C angular velocity {u) and average specific angular momentum 
on spherical shells (j). The enclosed mass at 10^° cm, where the jets were launched in our 
simulations, is 10.3 M©. 
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Fig. 2. — Time evolution of the density (bottom of each frame) and Lorentz factor (top of 
each frame) in the jet and its environs for Model 2A. The density is on a logarithmic scale, 
the Lorentz factor is on a linear scale, both color coded. Quantities are given 5, 10, 12, 20, 
40, and 70 s after the initiation of the jet at 0.1 x 10^^ cm. The third panel at 12 s is just 
as the jet is erupting from the star (radius 0.89 x 10''^''^ cm). By 70s, most of the internal 
energy has converted to kinetic and the Lorentz factor has almost reached its terminal value. 
Values of F at this time in the central, over-exposed jet are near 150. Note the explosion 
of the cocoon as the jet erupts and the ejection of material with terminal Lorentz factor, 
F ^ 10 — 20, at angles around 10 degrees off axis. 
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Fig. 3. — Time evolution of the density (bottom of each frame) and Lorentz factor (top of 
each frame) in the jet and its environs for Model 2B. The density is on a logarithmic scale, 
the Lorentz factor is on a linear scale, both color coded. Quantities are given 4, 8, 10, 18, 
40, and 70 s after the initiation of the jet at 0.1 x 10^^ cm. See also Fig. 2. Model 2B is a 
more energetic jet and reaches the surface in a shorter time. 
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Model 2C 
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Fig. 4. — Time evolution of the density (bottom of each frame) and Lorentz factor (top of 
each frame) in the jet and its environs for Model 2C. The density is on a logarithmic scale, 
the Lorentz factor is on a linear scale, both color coded. Quantities are given 8, 16, 18, 28, 
48, and 70 s after the initiation of the jet at 0.1 x 10^^ cm. See also Fig. 2. Model 2C is a 
less energetic jet and takes longer to reach the surface. 
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Model 2A 



55 F 




Angle 



Fig. 5. — Equivalent isotropic energy for Model 2A. As defined in the text, the equivalent 
energy to an isotropic explosion inferred by a viewer at angle 6 is plotted for various Lorentz 
factors. The line gives the energy contained in matter with V greater than the indicated 
value moving at a given angle. The top panel shows this quantity at 70 s; the bottom panel 
shows the estimated value much later, when all internal energy has converted to kinetic. The 
light gray line is a simple power-law fit, E-^^o = 1.5 x 10^^ x (^/2°)~^erg. 
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Fig. 6. — Equivalent isotropic energy for Model 2B. The light gray line is a simple power-law 
fit, Eiso = 4.5 X 10^^ X (^/2°)-3erg. See also Fig. 5. 
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Fig. 7. — Equivalent isotropic energy for Model 2C. The light gray line is a simple power-law 
fit, Eiso = 6.8 X 10^^ X ((9/2°)~^ erg. See also Fig. 5. The sohd angle between 8 and 10 degrees 
is 4 times that inside of 3 degrees. That is, the observers are more likely to be off-axis than 
on- axis. 
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2 4 6 8 10 12 14 

Angle 



Fig. 8. — Fraction of energy inside a certain angle to the total energy on the grid for (a) 
Model 2A, (b) Model 2B, and (c) Model 2C. Different lines are for material with different 
Lorentz factors. It is clearly shown that highly relativist ic material is confined to a small 
angle ^ 4°. For Model 2C, mildly relativistic (F > 5) material at larger angles {9 > 10) 
contains about 20% of its total in mildly relativistic energy. This fraction is smaller for 
Models 2 A and 2B. 
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(10^11 cm) 

Fig. 9. — Resolution study in two dimensions. The jets in Models 2B and 2T had identical 
parameters, but the calculation was carried out in cylindrical grids having different resolution 
(Table 4). Model 2T had the lower resolution. 
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Fig. 10. — A comparison of equivalent isotropic energy for Models 2B and 2T. The thin and 
dark lines show the equivalent energy for Model 2B at 10 s; the thick and light lines show 
the equivalent energy for Model 2T at 9.5 s. See also Tables 2, 4 and Fig. 9. 
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0.2 0.4 0.6 0.8 1.0 1.2 0.2 0.4 0.6 0.8 1.0 1.2 



z (10^11 cm) 



Fig. 11. — A comparison of two- and three-dimensional results at break out. Slices of the 
density distribution are shown along the polar axis. In the case of the three-dimensional 
models, the slice shown is along the x = plane. Model 2T and 3A had the same jet 
parameters and effective zoning (r and z in 2T, y and z in 3 A) and can thus be compared. 
Qualitatively the results at 9 s are similar, though Model 3A has evolved just a little bit 
faster. A dense plug of matter is visible at the head of the jet in both studies. Even though 
3A is a three-dimensional study, it retains the two-dimensional symmetry imposed by the 
jet's initial parameters. Model 3BS is like Model 3A, but with slightly asymmetric initial 
conditions, a 1% excess of energy in one half of the jet (see text). The asymmetry is such 
that a little more energy was deposited in the top half of the plane displayed than in the 
bottom half. The result is similar to Model 3A but now the strict two-dimensional symmetry 
is broken. The top of the jet looks quite different from the bottom. Model 3BL carries the 
symmetry breaking further with a 10% imbalance in energy at the base of the jet. The 
density of the "plug" is greatly diminished in these asymmetric jets. 
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Fig. 12. — Precessing jets - angle sensitivity study. Slices of the precessing jet models defined 
in Table 6 and § 4.1 corresponding to the x = plane are shown just after break out. The 
three jets had an inclination with respect to the radial of 3, 5, and 10 degrees, respectively, 
for models 3P3, 3P5, and 3P10 and the initial jet precessed on a cone with this half angle 
with a period of 2 s. For Model 3P3, the jet still emerges relatively intact (though one would 
want to follow the evolution further before concluding a GRB is still produced). The jets 
in Models 3P5 and especially 3P10 dissipate their energy before escaping and are unlikely 
to give terminal Lorentz factors as high as 200 (though they may still produce hard X-ray 
fiashes) . 



